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ABSTRACT 

We present bhlight, a numerical scheme for solving the equations of gen¬ 
eral relativistic radiation magnetohydrodynamics (GRRMHD) using a direct 
Monte Garlo solution of the frequency-dependent radiative transport equation, 
bhlight is designed to evolve black hole accretion flows at intermediate accretion 
rate, in the regime between the classical radiatively efficient disk and the radia- 
tively inefficient accretion flow (RIAF), in which global radiative effects play a 
sub-dominant but non-negligible role in disk dynamics. We describe the govern¬ 
ing equations, numerical method, idiosyncrasies of our implementation, and a 
suite of test and convergence results. We also describe example applications to 
radiative Bondi accretion and to a slowly accreting Kerr black hole in axisym- 
metry. 


1. Introduction 

Many of the brightest objects in the universe, including quasars and the lesser active 
galactic nuclei, stellar-mass black hole binaries, and gamma-ray bursts, are likely the results 
of black hole accretion driven at least in part by the magnetorotational instability (MRI; 
Balbus & Hawley 1991). The structure of the luminous plasma surrounding the black hole 
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remains uncertain (see the recent review of Begelman 2014), because it is difficult to resolve 
and because of physical complexity: relativistic gravity, turbulence in a magnetized plasma, 
and radiation transport all play some role in determining accretion flow structure. 

Nonetheless, accreting black holes may be partially classihed according to the ratio of 
their luminosity L to the Eddington luminosity, EEdd = 4:TtGMc/K es- 

For L > 10“^LEdd, radiation is dynamically important. Up to L ~ i^Edd, this regime 
can be modeled by the aligned thin a disk model of Shakura & Sunyaev (1973), in which 
the disk is geometrically thin and optically thick, and in which radiation pressure exceeds 
gas pressure at radii where most of the disk luminosity is produced. The radiative efficiency 
of the accretion flow, 77 ( 0 *) = L/(Mc^), is expected to be approximately constant and 
determined by the dimensionless black hole spin — 1 < a* < 1. It is common to describe the 
accretion rate M in units of an Eddington rate defined using a nominal efficiency rj = 0.1: 
m = rjMc^/L-Edd- For L > L^dd, the flow is expected to resemble the slim disk solution of 
Abramowicz et ah (1988), in which the flow becomes geometrically thick as a result of long 
radiation diffusion times. An obstacle to fully modeling the innermost, relativistic regions 
of flows with m > 10“^ is the need for an efficient relativistic radiation hydrodynamics 
scheme that can operate in both the optically thick (disk midplane) and optically thin (disk 
atmosphere, corona, funnel) regimes. 

For L UEdd, or m 1, accretion is likely to occur through a radiatively inefficient 
accretion flow (RIAF or ADAF; see the recent review by Yuan & Narayan 2014), in which 
the cooling time of a parcel of plasma is much longer than the time required for it to fall 
into the black hole. Radiation plays no role in determining the flow structure. RIAFs 
are believed to be geometrically thick, optically thin, collisionless plasmas that are at least 
partially supported by rotation. RIAFs are commonly modeled numerically using relativistic 
magnetohydrodynamic (MHD) codes, but it is unclear how well the fluid model describes 
the dynamics of the magnetized, collisionless plasma. It is also unclear how best to model 
the electrons, which are collisionally decoupled from the ions and determine the radiative 
properties of the plasma. However, local models, particularly numerical kinetic calculations, 
are beginning to constrain the electron distribution function in this regime (e.g. Kunz et ah 
2014, Riquelme et ah 2014, Sironi 2014). 

Between thin disks and RIAFs lies an intermediate regime in which radiation plays a 
modest role in the accretion flow; this configuration may be thought of as a RIAF perturbed 
by radiative effects. ADAF solutions evaluated at these accretion rates indicate a flow that 
is optically thin to Compton scattering (r ~ 10“^ — 10“^; Yuan et al. 2006), and optically 
thick only to synchrotron self-absorption at long wavelengths. As accretion rate increases 
the hrst non-negligible radiation-plasma interactions are expected to be Compton cooling 



and synchrotron cooling. For example, M87’s central black hole, an object of interest for the 
Event Horizon Telescope (Doeleman et ah 2009), is expected to reside in this intermediate 
regime (m < 6.3 x 10“®, based on a RIAF model; Kno et ah 2014), in Moscibrodzka et ah 
(2011) and Dexter et ah (2012). Snch systems exhibit nonlinear evolntion of conpled gas 
and radiation in strong gravity; predictive modeling is onr primary motivation for bhlight, 
a nnmerical scheme for general relativistic radiation magnetohydrodynamics. 

In the nonrelativistic and 0{v/c) regimes, many nnmerical methods have been developed 
to solve the radiation hydrodynamics eqnations (see the comprehensive review of Castor 
2004), inclnding flux-limited diffusion. Of particular relevance to black hole accretion flows 
is recent work on accretion in the near-Eddington regime using flux-limited diffusion (Hirose 
et al. 2009, 2014) and using the more accurate short characteristics method, in which specihc 
intensity is discretized in angle for each grid zone (Stone et al. 1992; Jiang et al. 2012, 2014a,b) 
and one obtains a full solution to the grey transfer equation. 

Close to the event horizon special and general relativistic effects can produce order unity 
variations in the intensity. These effects are particularly important for rapidly rotating black 
holes. Numerical schemes for solving the equations of general relativistic radiation MHD 
(GRRMHD) have only been developed in the last few years. All are frequency-integrated 
and use approximate closure schemes, including the Eddington approximation (Farris et al. 
2008, Zanotti et al. 2011, and Fragile et al. 2012) and low-order truncated moment closure 
(Shibata et al. 2011, Sqdowski et al. 2013, and McKinney et al. 2014). These schemes are 
formally accurate at high optical depth, but not for general flows at the modest optical 
depths relevant to black holes in the intermediate accretion rate regime. 

An alternative treatment of radiation, the Monte Carlo technique, has long been used 
for solving the full frequency-dependent transport equation without recourse to any closure 
model. Several radiation hydrodynamics schemes have recently been employed in astro¬ 
physics which couple a Monte Carlo representation of the radiation to a fluid model through 
interactions evaluated on a per-sample basis, yielding a Monte Carlo Radiation Hydrody¬ 
namics (MCRHD) scheme. This technique has received particular attention in the stellar 
physics community in Haworth & Harries (2012), Noebauer et al. (2012), Abdikamalov et al. 
(2012), Wollaeger et al. (2013), and Roth & Kasen (2014), which have variously investigated 
extensions such as implicit methods and interfacing the Monte Carlo representation with a 
continuum approximation in regions of large optical depth and/or large ratio of radiation 
to gas pressure, where the unadorned Monte Carlo technique fails. MCRHD schemes have 
also been implemented for studying star formation in Harries (2015), and have been used 
to model Compton cooling of accretion disks around black holes in flat space by Ghosh 
et al. (2011) and Garain et al. (2012). Monte Carlo techniques are particularly attractive 
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for GRRMHD because they are algorithmically simple, naturally incorporate frequency de¬ 
pendence (useful for treating Compton scattering) and the potentially complicated angular 
dependence expected in an optically thin regime, and are easily modihed to include special 
and general relativistic effects. 

In what follows we develop a scheme for GRRMHD called bhlight that is designed 
to model accretion flows with modest to low optical depth, bhlight couples two existing 
schemes: the GRMHD code harm^ (Gammie et al. 2003), and the Monte Carlo radiative 
transport scheme grmonty^ (Dolence et al. 2009). The paper is organized as follows: §2 
recounts the governing equations as they are solved in bhlight; §3 describes the numerical 
method; §4 demonstrates that bhlight converges on a set of test problems; §5 describes 
example applications to a radiating Bondi flow and an M87-hke disk model; §6 concludes. 


2. Basic Equations 

We adopt a physical model in which emission, absorption, and scattering of photons 
couple an ideal, magnetized fluid to the radiation held. We consider the huid and radiation 
sector in turn. The basic equations are identical to those integrated in the harm code (Gam¬ 
mie et al. 2003) and in the grmonty code (Dolence et al. 2009), but recounted here to dehne 
variables and expose physical assumptions. 


2.1. Fluid 

We assume particle number conservation, which in a coordinate basis is 

dt {y/^pou^) = -di {y/^pou") , ( 1 ) 

where po is the comoving frame rest mass density and is the huid four-velocity. 

Energy and momentum conservation for the coupled huid and radiation system are given 

by 

(T-t + = 0, (2) 

where is the magnetohydrodynamic stress-energy tensor, and is the radiation stress- 
energy tensor (not to be confused with the Ricci tensor). In a coordinate basis, Eq. 2 


^Freely available; http://rainman.astro.Illinois.edu/codelib/codes/harm/harm.tgz 
^Freely available; http://rainman.astro.illinois.edu/codelib/codes/grmonty/grmonty.tgz 
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becomes 

s, (V^n) = -a, (a=jrv) + + a=jG., 

where the radiation four-force density 






(3) 


(4) 


In the ideal MHD limit = 0 = electromagnetic held tensor), and one can show 

that 

n = (po + w + P + h‘^)u^u, + (P + ]f)g^ - (5) 

where P = huid pressure, u = huid internal energy density, and = b'^b^, with b^ the 

magnetic held four-vector, 

b^ = ( 6 ) 

and = — [nuK,\]/y^^ is the Levi-Civita tensor. 

Evidently b^u^ = 0, so b^ has only three degrees of freedom, expressed as P* = *p**, 
where 


Then 


*F^^ = = bV - b^u^. 


6 * = 

p* -i- pf 


p = 


w- 


The magnetic held evolution is determined by 

a, (v^B-) = d, [V^ (¥u‘ - bV)] 

subject to the no-monopoles constraint 

d, (v^P*) = 0. 

The equation of state is 

P = (7 — l)u. 


(7) 

( 8 ) 
(9) 

( 10 ) 

( 11 ) 

( 12 ) 


To summarize: the governing equations for the huid evolution are equations 1, 3, and 10, 
together with 11 and 12 . 
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2.2. Radiation 


The radiation field consists of photons with wave four-vector and momentum = 
hk^. The photons follow geodesics, with 






(13) 


and 

dk^^ 

^ ( 14 ) 

where T^^j, is the connection and A is an affine parameter along the geodesic. We assume 
that plasma dispersion effects are negligible, so photons travel on null geodesics, kfj_k^ = 0. 
The frequency of a photon in a frame with four-velocity is a; = —k^u^ {v = ci;/(27r)). 

In nonrelativistic radiative transfer one describes the radiation field with the specific in¬ 
tensity ly (here and throughout we ignore polarization), which is frame-dependent. However, 
oc fa where fa is the radiation distribution function 


where d^p = dpidp 2 dp 3 . Because diV, d^xp^^/^, and d^p/{p^y/^) are invariant, fu is also 
invariant. 


The evolution of fji is given by the Boltzmann equation, 

^ = CM (16) 

where A is an affine parameter along a photon trajectory (geodesic) and C[//j] accounts for 
interactions with matter: emission, absorption, and scattering of photons. The Liouville 
operator D/dA is a derivative along the photon trajectory in phase space. 

One can rewrite the Boltzmann equation as the radiative transfer equation. 


D 
dA V 



( 17 ) 


Here the extinction coefficient 


X.u — T CTi/, 


and the emission coefficient 

ffc = jv + v'AQ, 


(18) 

( 19 ) 
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where j^, is the fluid emissivity, is the scattering contribution to emissivity, is the 

scattering extinction coefficient, and is the absorption coefficient. Each of the quantities 
in parentheses in (17) is invariant. 

We neglect stimulated Compton scattering. The ratio of stimulated to spontaneous 
scattering is the photon occupation number in the scattered state. Models of highly sub- 
Eddington accretion onto supermassive black holes commonly feature: (1) relativistic elec¬ 
trons with ©e = kTe/{nieC^) > 1, corresponding to a mean amplification factor after Comp¬ 
ton scattering of ~ 160g; (2) a low frequency (millimeter or far-IR) peak in the spectrum 
at Vpk where the synchrotron absorption optical depth is 0{1). The energetically important 
single scattering events therefore produce scattered photons with z/^g ~ z/pfcl60g. Eor moder¬ 
ate accretion rates (i.e. scattering depth < 1 for the disk), the photon occupation number 
at Use is small, and so stimulated Compton scattering will be negligible. 

A consequence of our neglect of stimulated Compton scattering is that in a purely 
scattering medium the radiation field will approach a Wien (Boltzmann) distribution rather 
than a Bose-Einstein distribution. We verify this in Section 4.2. 

To summarize: the governing equations for the radiation are (17), (14), and (13), to¬ 
gether with appropriate expressions for the emission, scattering, and absorption coefficients. 


2.3. Radiation-Fluid Interactions 


In this section we adopt units such that c = 1 unless otherwise stated. It is apparent 
from Equations (3) and (17) that the fluid acts on the radiation through extinction and 
emission coefficients Xu ^md r]„. The radiation acts on the fluid through the four-force 
density We want to make these representations consistent. Begin with the manifestly 
covariant expression 

d^p 


Then 


= -V, 


RI,U ^ 


d^p 




P^p'^fR- 


:P^P^fR = -h 




-p- 


( 20 ) 


( 21 ) 


_ J v^p‘" dA 

The last equality follows from an expansion of D/dA and an integration by parts over mo¬ 
mentum space (Lindquist 1966). 

Using fu = and equations (17) and (21), 


G^ = — 


d^p 

v^P' 




( 22 ) 



where z/, Xu, lu, and 77 ^ are all evaluated in a frame with four-velocity u^. 

can be evaluated in an orthonormal tetrad frame comoving with the fluid 


^(a)’ ^(t) “ “ 


( 23 ) 


We will call this the “fluid frame.” In the fluid frame, 


G(a) = / dz/dfi {Xulu - Vu) ^(a), 


(24) 


where n(a) = p(^a)l{hv). Then 


G^^ = 

(a) 


(25) 


which is manifestly consistent with energy-momentum gains and losses by the radiation held. 


3. Numerical Method 

bhlight combines a second order hux-conservative ideal GRMHD integrator (Gammie 
et ah 2003) with a Monte Garlo scheme for radiation transport (Dolence et ah 2009) through 
radiation-huid interactions into a fully explicit GRRMHD scheme that is second order in 
space and hrst order in time for smooth hows. In this work we restrict ourselves to one- 
and two-dimensional hows, although the scheme can be trivially generalized to three spatial 
dimensions. 


3.1. Fluid Integration 

The huid integrator in bhlight is taken from harm, a conservative second order shock¬ 
capturing scheme on a two-dimensional mesh with an arbitrary spacetime metric. Here we 
give a brief summary of the method. Also, we adopt units such that c = 1, and for black 
holes we set GM = 1. 

The huid sector in bhlight updates a set of conserved variables U: 

U = ^(pou\T\,TiB^), (26) 

corresponding to the variables whose coordinate time derivatives are given in §2.1. These 
conserved variables are updated by huxes F; 


- B^v^) , 


¥=^g{p„u',T\,r,,B'W 


(27) 
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which in turn are calculated from the primitive variables P: 

p = (Po ,<•', B '), 


(28) 


where 


• 7 / 5 * 

= — 


(29) 


where n* = /vP is the fluid spatial 3-velocity, 7 = a/1 + giju^u^, a = is the lapse, 

and /5* = is the shift. Unlike n*, h* ranges over —cx) < h* < ex. In the Newtonian 
formulation all transformations between the nonrelativistic analogs of U, F, and P are 
analytic, but in the covariant formulation there is no general analytic form for P(U). 


The fluid update each timestep maps P” to its updated value P”+i by updating the 
conserved variables. Beginning with P”, the scheme calculates U** = U(P**) and F** = F(P**) 
via closed-form expressions, for F” after a reconstruction step that estimates P” at zone 
boundaries from values at zone centers. The update U"' —)■ over a timestep At is given 



where U represents the source terms such as those associated with the spacetime connection, 
values at n + 1/2 are estimated from a similar hrst-order step to and i,j here denote 

spatial indices in and respectively. This forms a second order, explicit timestepping 
scheme to t + At, and then P"-+i is found by numerically solving U(P"'’''^) = (see Noble 

et al. 2006 and Mignone & McKinney 2007). 

The fluxes F(P) are evaluated at zone faces using Local Lax Friedrichs fluxes. Prim¬ 
itive variables on either side of the zone face are determined through slope-limited linear 
reconstruction. We typically use the monotonized central limiter for reconstruction, but it 
is trivial to use higher order methods as well. 

Naively differencing the induction equation (10) will not preserve a numerical represen¬ 
tation of the no-monopoles constraint (11); the monopole density will undergo a random 
walk from zero with a step size determined by truncation error. Unless directly controlled, 
the monopole density can grow quickly and corrupt the solution. A variety of techniques 
for avoiding or removing magnetic monopoles exist; bhlight employs the flux-interpolated 
constrained transport (flux-CT) scheme introduced by Toth (2000). Although this intro¬ 
duces some additional diffusivity into the scheme, it is simple and effective. Details of the 
implementation are given in Gammie et al. (2003). 
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3.2. Radiation Transport 


bhlight uses nearly the same Monte Carlo implementation as grmonty, with a few 
important differences. The Monte Carlo samples are referred to here as superphotons. Each 
superphoton has a weight w (the number of photons carried by the superphoton), a momen- 
tnm = hkfj, and p^p^ = 0), and a position x^. 

The Monte Carlo representation of the the photon distribntion fnnction is 

fR,MC = S^{X^ - xl)S^{pj - Pj^k) (31) 

k 

where S^{x^ — x\) = 5{x^ — xl)S(x^ — xl)S(x^ — x^), etc. The snm is taken over all photon 
samples in the model, labeled by the index k, and Wk are the weights. Like /r, Jr^mc is 
invariant because Wk, S^{x^ ~ /{y/^p^), and S^{pj — Pj,k)V-^P^ (with pj covariant) are 
all invariant. 


Using eqnation (20), the stress-energy tensor is 




(32) 


This can be averaged over a three-volnme A^x = Ax^Ax'^Ax^ to obtain an estimate for 


Rf,U 


A^x 




E 


PkPk 

^Wk 


(33) 


,/^A^x^ PI 

where now the snm is taken only over photons within the three-volnme (zone) in qnestion 


3.2.1. Initializing the radiation field 

How should one initialize In bhlight’s target applications this question usnally 

does not arise because ffi relaxes rapidly to a quasi-equilibrium, so one can set = 0 in 
the initial conditions. In test problems, however, an accnrate initial Jr may be reqnired. In 
this case one wants to sample a set of photons in a single zone centered at Xc, that is, we 
want to sample A^xfR^Xc). 

One strategy is to sample jfi directly in a coordinate frame, using the invariance of 
Jr. For example, if jfi is thermal in the flnid frame, then the distribntion fnnction in any 
coordinate frame is where Bjj is the Planck function, and v = —u^k^/{2T[). 
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A second strategy, which we adopt, is to sample //? in the fluid frame (comoving 
indices are denoted by parentheses). Then we must take care: A^xffi{xc) is not invari¬ 
ant, because the volume element A^x is not invariant. But A^Xy/—gp^ is invariant, so 
[A^x/A^x') = y/—gp^/p^^\ where A^x' is a fluid frame volume element, and y/—g = 1 in the 
fluid frame. Then A^xf^lxc) = {A^x/A^x')A^x' Jr = {y/^p^/p^^^)A^x' Jr. This suggests 
that we can sample A^x'Jr in the fluid frame and then multiply the photon number diV by 
the corresponding —gp ^/to obtain a fair sample of A^x/r in the coordinate frame. 

This second strategy can be described more explicitly in terms of the Monte Carlo 
samples as follows. A list of photons in a single zone is obtained by taking 



fR^MC = '^'fXk 5^{Pj 

k 


Pj,k) 


(34) 


where the sum is over photon samples in a single zone. This is not invariant because S^{pj — 
Pj^k) is not invariant, so one must take care in sampling pj^k and Wk- Suppose we sample 
IrA^x' in the tetrad frame. This gives us a list of weights and momenta. We can transform 
back to the fluid frame using the invariance of \f-^p^d{pj — Pj^k), so that each WkS{pj — 
Pj^k) in the tetrad frame becomes Wk\/-^{p^/pd"^)5{pj — Pj^k) in the coordinate frame. We 
can therefore obtain a fair sample by adjusting the weights by a factor of ^/—in 
transforming from the fluid frame to the coordinate frame. 


3.2.2. Geodesic integration 

The position and wavevector of each superphoton is evolved individually by integrating 
Equations (13) and (14) numerically. Since evaluation of Christoffel symbols is costly, it is 
sensible to minimize the number of evaluations per timestep. 

We use the Verlet algorithm, a second order method that requires only one evaluation of 
the connection (number of evaluations is typically the order of the scheme). The algorithm 
as used in bhlight is identical to that used in grmonty and is described explicitly in Dolence 
et al. (2009). 

The Verlet method may be applied iteratively without re-evaluating the Christoffel sym¬ 
bols. For a fractional tolerance of 10“^ and the timesteps (corresponding to the AA) taken 
in bhlight, the scheme always converges. Although as of this writing we integrate all four 
components of k^, one could potentially integrate three components and use = 0 to 
evaluate the fourth, suppressing numerical errors and computational expense by a factor of 
4/3. 
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3.2.3. Units 


The radiation sector uses cgs units, except that photon wavevector components are 
measured in units of the electron rest mass energy. We therefore have to convert between 
units in the fluid sector and units in the radiation sector. For black hole spacetimes, this is 
done by choosing a cgs value for the fluid length unit and the fluid time unit, here 


C = 


GM 


(35) 


and 


T = 


GM 


(36) 


respectively. We also need a mass unit. Notice that the mass unit is not provided by the 
black hole mass in the test fluid {pC^ M) approximation used here. Instead we must scale 
the density, or equivalently the mass accretion rate, by choosing a cgs value for the mass 
unit M. Then, e.g., pcGS = PflvidM/C^. 


The components of the code photon wavevector are measured in units of nieC^. 
One might then be concerned about consistency between the transfer equation and the 
geodesic equation. The only condition for consistency is that the differential optical depth 
d-Ti, = which in turn requires that the ud\ = ds, i.e. that the units used in dehning 

u and dX be consistent and that the correct conversion be made from fluid sector units to cgs. 
In practice, then, we evaluate uKu in cgs units in the fluid frame and set dry = J\f{uK,u)d\, 
with Af = hC/{meC^). 


3.2.4. Superphoton Weighting 

The passive Monte Carlo code grmonty is designed to maximize the signal to noise in 
the hnal spectrum, which is measured in logarithmic intervals in frequency at spatial inhnity. 
The optimum allocation of weights would then place equal numbers of superphotons in each 
bin in log v. This requires an estimate of the hnal spectrum; grmonty estimates the hnal 
spectrum by integrating over the simulation volume, assuming the how is optically thin at 
all frequencies, neglecting gravitational redshift and Doppler shift, and setting the weights 
accordingly. 

In bhlight, by contrast, the weights should be designed to minimize errors in the 
dynamical evolution, i.e. in G^. The momentum and energy exchange associated with each 
radiation-huid interaction is proportional to wkn, where n is the huid frame frequency. This 
suggests that we should distribute energy uniformly among superphotons (as in Abbott & 
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Lucy 1985) to minimize the interaction noise, and thus set w oc 1/v. This is not generally 
possible because the four-velocity fluctuates across the simulation domain, but we will not 
do too badly if we ignore Doppler shift and gravitational redshift and thus set w ocl/v when 
sampling the emissivity. 

bhlight’s constant-energy-per-superphoton weighting scheme limits spectral resolution 
at low and high frequencies where the specihc energy density is small prior to scattering. 
These parts of the spectrum have little impact on the dynamical evolution, however, and 
higher quality spectra can be extracted in post-processing using grmonty. 


3.2.5. Emissivity 


At each timestep we sample the emissivity in the invariant four-volume y/—gAtA^x of 
each zone based on the fluid values at the half-step. It is easiest to sample the fluid emission 
in a comoving tetrad, where we have a simple expression for the emissivity. 


The emissivity is 




1 dE 
d^xdtdudil 


Using dE = hvdN where N is the number of photons, we can write 



h dN 
y/-^ d^xdtd log u 


(37) 


(38) 


Since dN = wdNg where Ng is the number of superphotons, we can then write for the number 
of superphotons produced per logarithmic interval in a single zone with volume A^x: 


dNg 

d\og V 


y/^AtA^X 


1 

hwiv) 



(39) 


In writing (39) we have made use of the invariance of y/—gd^xdt. Here w{u) oc 1/z/, and the 
constant of proportionality is set dynamically to keep the number of superphotons in the 
computational domain approximately constant. 


bhlight samples equation 39 between minimum and maximum frequencies Vrmn and 
^max, where the limits are set so that vji, is small outside this region in frequency space. It 
then uses a rejection scheme, sampling a uniform distribution in log < log v < log u^ax 
and a uniform distribution in 0 < r < (diVs/d log z/)jnax, where (diV 5 /dlogz/)max is the 
maximum of Equation 39. A sample is rejected when r > diV^/dlogz/. 

The angular distribution of photons is also sampled by rejection. The photon direction 
is given by where 6 is the angle between the magnetic held and photon direction in 
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the fluid frame and 0 is the corresponding azimuthal angle, bhlight samples a uniform 
distribution in 0 < cos 6 < 1, and a uniform distribntion in 0 < r < 1 . A sample is 
rejected if r > ju{0)/j^^rnax- It then samples a uniform distribntion in 0 < 0 < 2n. To 
ensure that the net force dne to emission in the fluid frame is zero to machine precision, 
photons are generated in pairs. Thus, a second photon is generated with the same frequency 
and cos 9' = —cos 9 and (j)' = <p + tt. In the fluid frame, = u, = a; sin 6 ^ cos 0, 
k^y) = casing sin 0, and k^^'^ = ujcos9^ where is parallel to the magnetic held. Once 
we have a snperphoton sample in the comoving frame it is transformed to the coordinate 
frame using a pre-constrncted orthonormal tetrad The snperphoton is set to the 
zone center to avoid additional orthonormal tetrad construction.^ 

Sampling is a snbdominant compntational expense in bhlight, so althongh one could 
develop more efficient sampling schemes, a simple rejection scheme is adequate. 

Four-momentum is locally conserved and so snperphoton emission implies a back-reaction 
on the emitting fluid. A pair of emitted superphotons with wavevectors /cj*, k^ correspond 
to a change in fonr-momentum (in flnid code units): 

= Vw {k't + k^) , (40) 


where V = rrie/Ai, which in tnrn specifies the contribution to the four-force density AG^: 


AG, 




(41) 


where all geometric quantities are evaluated at zone centers. Because photons are emitted 
in pairs, the spatial components of k^ cancel in the fluid frame, and Sp, oc u,. 


3.2.6. Absorption 

grmonty treats absorption deterministically by continuonsly decrementing w along a 
ray. A similar deterministic procedure has been shown to suppress noise in Monte Carlo 
radiation hydro schemes (e.g. Noebaner et al. 2012) in flat space. However, formulating snch 
a scheme in general relativity, where photons move along geodesics, is more complicated 
because the photons follow cnrved trajectories through each zone. 


^Because photons are created at zone centers, our scheme will fail when individual zones become optically 
thick. Should this become a problem the scheme can be modified so that new superphotons are distributed 
within a zone. 
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In bhlight we treat absorption probabilistically. While stepping a superphoton by AA 
along a geodesic, the incremental optical depth to absorption Ar^ = A/'/tj^ab^z/AA. Here 
T^i^u,abs is the invariant absorption coefficient, evaluated in the fluid frame and interpolated 
to x^. An absorption occurs if 

Ara>-logra, (42) 

where 0 < < 1 is sampled uniformly; the absorption occurs at Axa = — logra. To process 

the event we push the superphoton back, A —)■ A + AA(logra/Ara), put the superphoton 
four-momentum into the fluid at that location, and annihilate the superphoton. 


The four-momentum change in the fluid Ap^ due to absorption of a superphoton with 
wavevector is 

Ap^ = VwkA (43) 

This can be expressed as a contribution to the radiation four-force density AG^: 


A — 

^ ^A^xAf 


(44) 


where y/—g is evaluated at the zone center. 


3.2.7. Scattering 


We treat scattering probabilistically in bhlight, as in grmonty. Scattering is similar to 
absorption, i.e. scattering occurs when 

Ats>-{\ ogrs)/bs, (45) 

where 0 < r* < 1 is sampled uniformly, Ar^ = Mk^^s^AX, and is the invariant scattering 
opacity. Because scattering events are rare but energetically important we have introduced 
a bias parameter 6^ > 1 to enhance the probability of sampling scattering events. To process 
the event, we push the photon back along the geodesic from A-|-AA to A-|-AA logrs/(6sArs). 
To preserve photon number a scattered superphoton is created with weight Ws = w/hs and 
the original superphoton has weight set to w' = w — Wg. 


In general, a superphoton is subject to both absorption and scattering simultaneously. 
In a deterministic treatment, the code must dynamically choose which process, if any, to 
apply to the superphoton. To handle this in an unbiased manner, for each photon we, 
assuming that at least one of the inequalities Eqns. 42 and 45 has been satisfied, choose 
which interaction to process according to a similar weighted sampling. That is, for 


- log^g 
AXa 


- log Ts 

bsAxs 


< 


( 46 ) 
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the absorption interaction is chosen; else, scattering is chosen. With this method, large 
optical depth or bias in one interaction will not serve to decrease the physical effect of the 
other interaction, although it will increase the number of superphotons required to resolve 
both interactions simultaneously. 


How should we set 6^? Most of our models have 1, so only superphotons 

would produce scattering events if bg = 1, and the energy per superphoton would increase 
by the mean amplihcation factor H 1 + 40e + 160g.^ This suggests that we should set 
bs ^ A to maintain constant energy per superphoton. There are two failure modes to be 
avoided, however. First, if 6^ ^ f/Zs then each superphoton will scatter more than once 
and the number of superphotons on the grid will grow exponentially. Second, if Awhu is 
larger than the total internal energy in the zone where the scattering occurs then the zone 
energy will be negative after scattering (more on this in the next subsection). Together, this 
suggests that we set 


bs = A MAX C, 


. (47) 

u^A^x) ^ ’ 

Here C depends only on t and is dynamically adjusted to control the number of scattered 
superphotons in the simulation. The requirement bgTg ~ 1 is equivalent to each superphoton 
scattering once between emission and escape through the boundary (neglecting absorption). 
Over a timestep At, one can estimate the number of photons which escape through the 
boundaries of a domain with linear dimension L as NcAt/L, where N is the desired total 
number of superphotons (which sets the weights for emission as described previously). To 
enforce the requirement that each superphoton scatter once, C is calculated dynamically as 
the ratio of this estimate to the real number of scattering events per timestep, averaged over 
some timescale. 


Each scattered superphoton is generated from an incident superphoton wavevector 
as follows. The four-momentum of the scattering electron is sampled from a thermal 
(Maxwell-Jiittner) distribution according to the procedure described in Canheld et al. (1987). 
The scattered superphoton wavevector is sampled from the Klein-Nishina differential 
scattering cross section in the rest frame of the scattering electron and boosted to the fluid 
frame and then transformed to the coordinate frame. It is then assigned a weight and entered 
in the list of active superphotons. 

Each scattering event generates a change in fluid four-momentum, 

Ap^ = (F - A;,^) (48) 


^This approximate expression overestimates A by 16% at ©e « 1/2. A better estimate, which underesti¬ 
mates A by 4% at ©e « 0.02, is A — 1 « 4©e — 2©e^^ -I-16©^. 
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and a corresponding contribntion to the flnid throngh the fonr-force density 


AG,= 




where y/—g is evalnated at the zone center. 


(49) 


3.3. Radiation Force in Fluid Evolution 

The radiation force is treated in an operator-split fashion. The flnid integrator initially 
npdates the conserved variables U from step n to n -|- 1 over the entire timestep At without 
radiation, i.e. it performs U” —)■ as described in Section 3.1. This fluid integration 

generates half-step fluid primitive variables P"+i/2; these values are sent to the radiation 
sector and used to evaluate the total radiation four-force density for each zone. 

The fluid integrator then updates the fluid variables with the radiation interaction t 

U«+i tjy considering only the radiation contribution to the conserved energy and momentum 
variables: 

+ A«v/=5G., (50) 

jjn+i gj-Q ^]^g gnal conserved fluid variables at the (n -|- 1)*'^ step. The evolution is 
therefore hrst order in time. 

The evolution is explicit and the radiation and fluid share a common timestep At, which 
we set to the minimum grid zone light crossing time ~ Axjc, where Ax is a characteristic 
zone lengthscale. As is well known, the radiation source terms are stiff when the timescale for 
exchange of energy-momentum between the fluid and radiation is smaller than a timestep. 
The cooling time Tcooi = u/A, where A = cooling rate per unit volume, = u^G^ ~ Uradc/^mfp, 
where Xmfp is a suitably frequency-averaged absorption mean free path and Urad = 
is the radiation energy density in the fluid frame (one can perform a similar estimate for 
Compton cooling). Thus the source term is stiff if u/A < Ax/c or {urad/u){Ax/Xmfp) > 1, 
or when the optical depth across a zone exceeds u/urad- 

For our scheme we must also consider robustness in the presence of Monte Carlo noise. 
Even if Tcooi/A t > 1 the cooling rate may fluctuate upward so that a zone loses all its thermal 
energy in a single timestep. This can happen if m/(A^)^A < Ax/c. This will differ from 
the usual stiffness condition only when the number of absorption events per timestep in 
a zone is small compared to one. The condition for robustness against this failure mode, 
“supercooling,” where a single photon causes the zone to lose all its internal energy, is that 
whv < uA^x (where we have left out geometric factors). 
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Where, then, will bhlight fail? The radiation force source terms are stiff when the 
optical depth across a single zone exceeds u/urad- For black hole accretion applications we 
expect this only for models in the high accretion rate regime, M > 10~^MEdd, although 
the precise condition will depend on details of the evolution and the numerical setup. This 
problem could be remedied by using an implicit update, but Monte Carlo is probably not 
the optimal method for studying this regime anyway. The supercooling problem is more 
relevant for our target application to intermediate accretion rate black holes, and arises if 
the internal energy content of a zone is small compared to the typical superphoton energy. 
This can occur in low density regions over the poles of the black hole, but the fluid evolution is 
inaccurate there in any case (because the truncation error in internal energy is dominated by 
the magnetic held evolution, to which it is coupled via the total energy density), and negative 
internal energies are dealt with by harm’s hoor routines, resulting in a small nonconservation 
of energy. 


3.4. Parallelization 

bhlight is a hybrid MPI/OpenMP scheme in which a single node handles the huid inte¬ 
gration, multiple nodes evolve the radiation, and a single additional node acts as gatekeeper 
between the huid and radiation sectors. During each timestep, the only exchanges are an 
array of radiation four-force densities to the huid sector and an array of huid variables to the 
radiation sector via the gatekeeper node. The gatekeeper node distributes the huid variables 
to all radiation nodes, and reduces the four-force density contributions from each radiation 
node. 

We evolve radiation on each node independently from other radiation nodes. After 
globally scaling the emission weights and scattering bias to yield approximately the desired 
number of superphotons at saturation, the code samples emission events on each radiation 
node, each of which has access to the entire array of huid variables and maintains its own set 
of superphoton samples. Emission, absorption, and scattering events generate a four-force 
density contribution. At the end of every timestep, these contributions are reduced by the 
primary radiation node over MPI. Each radiation node is individually parallelized under 
OpenMP, further dividing the superphoton calculations across individual compute cores. 
We parallelize the main compute loops for the huid sector with OpenMP, which enables 
completion in a reasonable clock time for an axisymmetric calculation. 



- 19 - 

3.5. Implementation Details 

In attempting to describe our numerical implementation in a coherent narrative we have 
omitted certain secondary topics, which we now collect here. 

• The radiation sector in bhlight makes extensive use of random numbers. We use the 
Mersenne Twister algorithm from the GNU Scientihc Library, with a different random 
seed for each MPI node and each OpenMP thread. 

• In axisymmetric disk calculations, we implement a form of static mesh rehnement by 

using modihed Kerr-Schild (MKS) coordinates {t,x^,x^}. and are related to 
the Kerr-Schild r, 6^ by r = exp -|- ro and 6 = -|- ((1 — hs)/2) sin(27ra;^), where 

ro G [ 0 , cx)) and hg G ( 0 , 1 ] are free parameters. 

• Truncation error in the geodesic integration causes to drift off the lightcone (this is a 
consequence of our decision to integrate all four components). We destroy superphotons 
with negative frequency in the fluid frame; for torus runs as in Section 5.2, we hud ~ 
1.1 X 10“® destructions per geodesic update. This problem does not occur in Cartesian 
coordinates in Minkowski space. 

• Scattered superphotons may scatter any number of additional times during the same 
timestep, provided sufficient optical depth to do so. 

• bhlight does not conserve momentum and energy to machine precision because of 
truncation error in the geodesic integrator. However, for the integrator tolerance and 
typical superphoton resolutions this is not signihcant (for torus runs as in Section 5.2). 
At some time, the average fractional error in energy relative to the initial energy at 
emission is ~ few x 10 “^; if this were to become a leading source of error, increasing 
the integrator tolerance is not a signihcant expense. 


4. Test Suite 

We have developed a suite of test problems for bhlight. Since the huid and radiation 
sectors of bhlight use well tested codes, we focus on problems with coupling between the 
two sectors. Good test problems are hard to hud, since there are few known exact solutions to 
the equations of radiation MHD in either Newtonian or relativistic contexts. We substitute 
approximate solutions to the full equations, such as some of the shocks we consider below. 
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We do not consider pure transport tests that are trivially satisfied by a Monte Carlo scheme, 
such as shadow tests, expanding pulses, and dynamic diffusion®. 


4.1. Optically Thin Cooling 


We consider the temperature evolution of an optically thin, radiating, stationary, ideal, 
and homogeneous gas initially at temperature Tq. The gas obeys a 7 -law equation of state 
i.e. p = (7 — l)u. The density and velocity of the gas are fixed; only the temperature is 
allowed to evolve. The bremsstrahlung-like emissivity is 

ji. = exp (—/izz/ZcbT), (51) 


where N = 5.4 x 10 cm® s ^ Sr ® Hz ^ is a constant. The associated cooling rate A 
is given by 


A = —— = [ dudilji, = dvr 
dt 


ksn^N 


j^l/2 


(52) 


which implies a temperature evolution 


m = To 1 1 - i 
*/ 


(53) 


1 /2 

valid from ti = 0 to tf = HTq /((y — l) 7 rA^n), the time at which the temperature of the 
fluid reaches zero. 

For this realization we choose 7 = 5/3, tf = 10® s, and Tq = 10® K. Figure 1 shows the 
resulting numerical evolution plotted against the analytic solution. Convergence in the 
norm, 

~ (54) 


(/*) — ^ ^ |/num,i freferenced 


_^ /2 

is expected to scale as Ns , this can be seen in Figure 2, which is evaluated near t = tf. 


4.2. Compton Cooling 

Consider a closed, one-zone model in which Compton scattering is the only permitted 
interaction between an ideal, 7 = 5/3 gas initially at temperature Tg^i and a swarm of 


^Such tests can be performed with the freely available grmonty code. 



21 


Fig. 1.— Optically thin cooling of one static fluid zone. Approximately 5 x 10® superphotons 
were created. 



photons all with initial frequency u = uq. Fluid motion is suppressed; only the internal 
energy is allowed to evolve. The number of photons is conserved, and in thermal equilibrium 
Tgj = Trj = Tf, the radiation approaches the Wien distribution, 

exp{-hu/{kBTf)), (55) 

where is the number density of photons. This tests the scattering kernel and Compton 
heating and cooling of the fluid. 

We set the electron number density n = 2.5 x 10^^ cm“®, Tg^i = 5 x 10^ K, uq = 
3 X 10^® Hz, and = 2.38 x 10^® cm“®. The characteristic (Compton) relaxation time is 
the photon mean free time (ngcrrc)"^ divided by the fractional energy change per scattering, 
~ kT/{meC^), yielding a Comptonization timescale ~ 0.02 s. We can predict the hnal state 
using: (1) thermal equilibrium; (2) conservation of photon number; (3) that the hnal photon 
distribution is Boltzmann; (4) conservation of total energy. This yields Tf = 5.19 x 10® 
K. Figure 3 shows huid and gas temperatures equilibrating at the correct temperature on 
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approximately the estimated timescale, and u^, = d-E/(d^a;dlogz/) plotted against the an¬ 
ticipated distribntion along with associated residnals. Note that we dehne T^. = ErKSn^ks) 
[Er is the radiation energy density), which assumes a Wien distribution, and so this value 
is strictly valid only at late time. 

4.3. Linearized Transfer and Energy Eqnations 

As a test of the full transfer equation in ID with gray absorption, consider a sinusoidal 
temperature perturbation in a static gas. Mihalas & Mihalas (1984) show that in this case the 
full transfer equation plus gas energy equation admit damped solutions. The eigenmode has 
perturbations oc exp{i{kx — t/tRn)) for wavenumber k and decay time tRR. The dispersion 
relation is 

(?))]■'• <-) 

where Tq is the mean temperature, p is the material density, ao is the frequency-integrated 
extinction coefficient, and c„ is the specific heat capacity. We simulate this problem in 
bhlight by initializing one wavelength in a ID box in local radiative equilibrium with 64 
grid zones and periodic boundary conditions. The amplitude of the initial perturbation is 
O.OSTq. We evolve this system for a variety of optical depths per wavelength r by varying 
ao. We obtain a decay time from the amplitude of the best £t sinusoid at t = tRR. We 
find good agreement with Equation 56 in both optically thin and optically thick regimes, as 
shown in Figure 4. 

We also examine convergence for r = 1 (with 100 grid zones); this result is shown in 

1 /2 

Figure 5. Evidently the errors scale as W' , as expected. 

4.4. Relativistic Radiation MHD Linear Modes 

We now also consider linear modes of the full equations of one-dimensional radiation 
magnetohydrodynamics; that is, we now include momentum exchange and magnetic fields. 
Acquiring even a linear solution to these equations with full transport is challenging, and so 
we resort to the approximate, relativistic Eddington closure scheme of Farris et al. (2008) 
with gray opacity k from which to extract plane-wave solutions P. Details of our calculation 
for a perturbation hP for 6 ~ exp (cat — ikx) about a thermal equilibrium Pq are given in 
Appendix C. With a magnetic held R* = (Rq, Bqi 0) we follow the nonrelativistic treatment 
of Jiang et al. (2012) by confining variation to a plane (suppressing Alfven waves): P = 
(p,n, u\u^, R^ R2, R, E\ uo + 6 u, 6u\ Rq, Rq + JR^ Ro + JR^ JR^). 
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Fig. 3.— Evolution of the Comptonization problem. The top panel shows radiation and 
gas temperatures approaching the analytic final temperature. The middle panel shows the 
initial and hnal Uy for the radiation, along with the analytic result for the hnal state. The 
bottom panel shows the residuals for the hnal spectrum; the numerical spectrum is apparently 
unbiased even at frequencies with low sampling resolution (shaded regions). 
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Fig. 4.— Dispersion relation for eigenmode of the transfer and energy equations as a function 
of optical depth per wavelength. Solid line shows analytic expectation, while points show 
bhlight results. At this resolution, the fractional error is ~ 10“^. 
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bhlight is not designed to evolve perturbed equilibria; we focus on cases which accomo¬ 
date both bhlight’s numerical limitations as well as the discrepancy between the Eddington 
closure and bhlight’s full transport (i.e. we consider only many optical depths per wave¬ 
length). We study convergence of two specihc cases with signihcant radiation pressure: a 
nonrelativistic radiation-modified slow MHD mode, and a relativistic radiation-modified fast 
MHD mode. In this section, c = = 1- For all calculations, we use a box of length L = 1 

with 128 grid zones, evolve to hnal time tf, set the wavenumber of the perturbation k = 27r, 
and normalize the to be < 1% of the Pq for all P.® The ratio of radiation to gas pressure 
Pr = and the optical depth per wavelength r = upL. 


4.4-1- Radiation-Modified Slow Mode 

For a magnetized fluid in the presence of radiation, the MHD modes are damped in a 
similar fashion to the radiation hydrodynamic case. We first consider the radiation-modihed 
slow mode solution. We set 7 = 5/3, p = 1, u = 0.01, Bq = a / 5 / 6 , fir = 1, and r = 20 
and evolve the initial conditions in bhlight to tf = 2.5, nearly half an e-folding time. The 
eigenmode is given in Table 1. Expected convergence at in the average number of extant 
superphotons Ng for Monte Carlo-dominated error is shown in Figure 6. 


4-4-2. Radiation-Modified Fast Mode 

We now consider the radiation-modihed fast mode solution, for a relativistic equilbrium. 
We set 7 = 4/3, p = 1, M = 10, Bq = a / 5 / 6 , fir = 1, and r = 20 and evolve the initial condi¬ 
tions in bhlight to tf = 1.7, approximately a wave period. Note that radiation damping is 
not signihcant during this time. The eigenmode is given in Table 2. Expected convergence 
at t/ in the average number of extant superphotons Ng for Monte Carlo-dominated error is 
shown in Figure 7. 


4.5. Su-Olson Problem 

Su & Olson (1996) found a solution in terms of integrals to the coupled energy balance 
and radiative transport equations in the dihusion approximation for a semi-inhnite slab of 


®The SageMath notebook used to evaluate these modes may be accessed via SageMathCloud at 
http://bit.ly/lCCi82y 



Fig. 6.— Convergence for the radiation-modified slow mode. 
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Table 1. Radiation-Modified Slow Mode 


u -0.155954250795 + 0.506371984839? 

~6p 0.992522043854 

5u 0.0115437955397 + 0.00253571930238? 

5u^ -0.0799889439467 - 0.024635280384i 

(5??2 -0.0804556011602 - 0.0252291311891? 

5B^ -0.00672014035309 - 0.00465692766557? 

6E 0.0129233747759 + 0.0201394332108? 

5F^ 0.00205715652365 - 0.00136719504591? 

6F^ -3.27455464963 x 10-^ + 5.02957595074 x lO"®? 
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Table 2. Radiation-Modified Fast Mode 


a; -0.000695187092855 - 3.6761842859* 

~Tp 0.0528655837266 - 0.000203781373769* 

5u 0.704834313006 - 0.00203965222393* 

5u^ 0.0309307265333 - 0.000125078175647* 

5u^ -0.00191512771962 + 0.000183786243559* 

5B‘^ 0.0512475714061 - 0.000472212042911* 

5E 0.704838422737 

5F^ -4.41864621469 x lO"® + 0.00198608271502i 

5F‘^ 0.000397994468532 - 0.00462049532989* 


Fig. 7.— Convergence for the radiation-modified fast mode. 
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static, initially cold fluid (with heat capacity c„ = with gray absorption coefficient a, 

and a Marshak (isotropic incident radiation) condition at the left boundary with incident 
flux F. The solution is given in terms of dimensionless gas and radiation energy densities u 
and V, respectively, versus the dimensionless spatial coordinate x = y/Saz and dimensionless 
time coordinate t = Aancat 'where an is the radiation constant, u and v are defined as: 



where E is the radiation energy density. 

In replicating this solution with bhlight, we adopt parameters such that the solution 
remains optically thick, without being so optically thick across a grid zone that our Monte 
Carlo transport scheme fails. We use 1024 grid zones on a; = [0, SOy^] and evolve the system 
to t = 500 (when the solution is approximately in equilibrium). Results are shown in Figure 

8 . 


4.6. Radiative Shocks 

We now turn to dynamical tests: here, ID radiative shocks. We will use nearly the 
same suite of relativistic, radiative shocks considered by Farris et al. (2008) in testing their 
GRRMHD code, which was based on a nonequilibrium Eddington closure. Because we solve 
the full transfer equation, we expect disagreement on scales of order the photon mean free 
path. 

The Farris et al. tests assume a grey opacity k and are set in Minkowski space. The so¬ 
lutions are described by po, the x component of the fluid four velocity the gas pressure P, 
and the comoving-frame radiation energy density E and x-component of the radiation flux 
four-vector F^. The latter vanishes far from the shock. We consider only the shock-frame 
(unboosted) version of the tests, initialized as shock tubes. All tests are purely hydrody¬ 
namic: the magnetic held plays no role. 

Our tests are identical to those given in Farris et al. (2008) except that we modify case 
(4), which has radiation pressure a factor of 10 larger than gas pressure upstream from the 
shock. Due to this large radiation pressure, bhlight cannot integrate this case stably with 
the numerical resources available to us. Instead, we set the upstream radiation pressure 
equal to the gas pressure, and call this case (4a). The shock parameters are listed in Table 3. 
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Fig. 8.— Numerical gas and radiation energy densities versus analytic gas and radiation en¬ 
ergy densities. The bhlight calculation and Su-Olson results show excellent correspondence 
at this late time. 
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Units are such that c = 1; (equivalently, h) is determined by enforcing thermal equilibrium 
E = aR^P/P qY far from the shock. 

Case 1 is a nonrelativistic strong shock with gas pressure much greater than radia¬ 
tion pressure, and consequently the fluid variable prohles resemble a nonrelativistic shock, 
bhlight output and the analytic solution are shown in Figure 9. Correspondence is good 
except for small deviations in the radiation variables near the shock interface, as expected 
for our full transfer solution. 

Case 2 is a mildly relativistic shock with somewhat larger radiation pressure than in 
Case 1. bhlight output and the analytic solution are shown in Figure 10. The prohles of E 
and E^ show qualitative differences between the bhlight result and the analytic solution: 
a discontinuity in the analytic solution does not appear for the case of full transfer. This is 
unsurprising given the approximate nature of the Farris solutions. Note that for this case 
the approximate Eddington solution contains an unphysical discontinuity in the coordinate 
frame radiation energy density. 

Case 3 is a highly relativistic shock with dynamically important radiation held, bhlight out¬ 
put and the analytic solution are shown in Figure 11. We hnd signihcant diherences within a 
few photon mean free paths of the shock, particularly for the radiation hux. Figure 12 shows 
the expected self-convergence in the radiation variables. A similar self-convergence 

trend also appears in the huid variables until grid resolution becomes the dominant source 
of error. 

Case 4a is a modestly relativistic wave with upstream radiation and gas pressure nearly 
equal, bhlight output and the analytic solution are shown in Figure 13. We hnd good 
agreement in all variables, despite the relatively strong radiation held. Note also that even 
with a large number of samples it is difficult to suppress noise in E^ when it is much smaller 
than E. 


4.7. Black Hole Atmosphere 

Next we turn to a general relativistic equilibrium test. Consider a Schwarzschild black 
hole surrounded by a static atmosphere. The atmosphere is bounded by static, concentric 
spherical shells at r* > 2GM/c^ and rp > r'i, and is in radiative equilibrium. The shells are 
rehecting boundaries and exchange no heat with the atmosphere, which has a grey opacity 
K and adiabatic index 7 = 5/3. 

In the Newtonian limit radiative conduction would drive the atmosphere toward T = 
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Table 3. Parameters for Farris shocks 

Case 7 n Left State Right State 

1 5/3 0.4 po = 1-0 Po = 2.4 

P = 3.0 X 10-5 P = 1.61 X 10“^ 
= 0.015 = 6.25 X 10-5 

P = 1.0x10-5 P = 2.51x10-^ 
5/3 0.2 Po = 1.0 Po = 3.11 

P = 4.0 X 10-5 P = 0.04512 
= 0.25 = 0.0804 

P = 2.0 X 10-5 E = 3.46 X 10-5 
3 2 0.3 Po = 1.0 Po = 8.0 

p = 60.0 P = 2.34 X 105 

= 10.0 = 1.25 

E = 2.0 P = 1.14x105 

~4a 5/3 0.4 po = 1.0 po = 1.165 

p = 0.1 P = 0.1233 

= 0.5 = 0.4292 

P = 0.3 P = 0.3763 
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Fig. 9.— Gas and radiation variables for radiative shock Case 1. The result from bhlight is 
shown as a red solid line, and the analytic solution is shown as a dashed line. 
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Fig. 10.— Gas and radiation variables for radiative shock Case 2. The result from bhlight is 
shown as a red solid line, and the analytic solution is shown as a dashed line. 
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Fig. 11.— Gas and radiation variables for radiative shock Case 3. The result from bhlight is 
shown as a red solid line, and the analytic solution is shown as a dashed line. 
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Fig. 12.— Self-convergence of all variables in radiative shock Case 3. Ng is directly propor¬ 
tional to the number of extant photons in the simulation. Dashed lines correspond to the 
Ns ' trend expected for Monte Carlo integration in the absence of resolution errors in the 
hydrodynamics solver. 
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Fig. 13.— Gas and radiation variables for radiative shock Case 4a. The resnlt from 
bhlight is shown as a red solid line, and the analytic solntion is shown as a dashed line. 
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const. In a relativistic atmosphere the gravitational held causes the atmosphere to come into 
a different equilibrium in which redshifted temperature, T^o = T^y—gQQ, is constant. 

In the Eddington approximation the atmospheric structure is determined by the con¬ 
ditions Too = const, and mechanical equilibrium, = 0. Note that any radiation source 
terms vanish in thermal equilibrium; in fact, the Eddington approximation solution turns 
out to be an exact solution to the transfer equation, as we show in Appendix B. 

Once the inner and outer radii are specihed, there are three dimensionless parameters 
that describe the solution (although their interpretation is purely Newtonian): the ratio of 
the atmospheric scale height at the inner boundary to the local radius hi = ’’^ksTiri/{iinipGM) 
(where fi is the mean molecular weight), the ratio of radiation pressure to gas pressure at the 
inner boundary I3r = fimpauT^/(SpikB), and the characteristic optical depth r = Kpi{ro — ri). 

We set Vi = 3GM/c^, To = 20GM/c^. We set h ^ 2.66, I3r ~ 0.23, and r ^ 5.0. We use 
a ID domain with 128 zones; the errors are dominated by Monte Carlo noise. The solution is 
shown in Figure 14. We also hnd the comoving spectrum at hducial radius r r* to match 
our expectation for thermalized radiation (with free vertical scale), as shown in Figure 15. 
Evidently the redshifted temperature is indeed constant. 
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Fig. 14.— Rest-mass density and both gas and radiation temperatures for the static atmo¬ 
sphere test at t = 150M. Dashed lines represent the analytic solutions. 
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Fig. 15.— Spectrum of superphotons at r r* for the radiating atmosphere test, showing 
good correspondence to the expected Planck spectrum. 
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5. Applications 
5.1. Radiating Bondi Accretion 

Here we report a preliminary investigation of spherically symmetric accretion onto a 
Schwarzschild black hole with radiative coupling. The fully dynamical problem has been 
studied analytically (although not with a full transport solution, frequency dependence, 
or magnetic helds) by Vitello (1984), Park (1990), and Nobili et ah (1991). Frequency- 
integrated numerical studies have also been performed in Fragile et al. (2012), Roedig et 
al. (2012), and McKinney et al. (2014) over many GMjc^. Here we obtain a full transport 
solution with frequency dependence and magnetic fields. We include synchrotron emission, 
synchrotron absorption, and Compton scattering, with associated heating and cooling. 

We set M = 6.6 x 10 ®Mq and describe the spacetime with modihed Kerr-Schild (MKS) 
coordinates (Section 3.5). All simulations are performed in ID, with r G [1.5M, 50M] for 64 
zones. As initial conditions for the fluid, we adopt the nonradiative Bondi solution of Hawley 
et al. (1984) (which we hold constant at the outer boundary), except we set 7 = 13/9 and 
place the sonic point at r = 200M. We vary the accretion rate by varying the density (or 
equivalently the mass unit Ai) of the flow. 

The magnetic field is initialized as a radial, monopolar field, which has no effect on the 
fluid motion. We set = a/r^, where a is chosen such that = 2P/lP‘ ^ 130 at r = 2M in 
the initial (GRMHD) conditions. No radiation is present initially. The radiation is allowed 
to flow freely out of the computational domain at the radial boundaries, with no inflow of 
radiation. 

The luminosity L{r), evaluated once the flow has settled and become almost time- 
independent, is 

L{r) = j (59) 

We average L{r) over radial shells between r = dOGM/c^ and r = 50GM/c^ to obtain an 
average luminosity L. This is equivalent to time-averaging at a single radius. 

We perform hve calculations at different accretion rates. We evolve each system until 
the luminosity becomes stable. At high accretion rate this can take as long as 400GM/c^. 
The characteristic number density for each simulation, mass accretion rate 

M = — J dx^dx^pu^ (60) 

(evaluated at the event horizon), and L are given in Table 4. The resulting prohles of 
these hve cases, along with that of Case 5 with Compton scattering disabled, are shown in 



43 


Figure 16. The lack of cooling in the pure synchrotron case indicates that Compton cooling 
dominates over synchrotron cooling for the highest accretion rate model. The relationship 
between luminosity and accretion rate when Compton scattering is active is shown in Figure 
17. 

We have checked self-convergence of a solution with n = 1.51 x 10^'^ cm“^ (between 
Cases 4 and 5 in Table 4) in steady state at t = 200M. The expected convergence behavior 
for Monte Carlo-dominated error is shown in Figure 18. 

Our models show that Compton cooling is important for Bondi accretion near the Ed¬ 
dington rate, and that - for our assumed held conhguration - synchrotron cooling is compar¬ 
atively unimportant. Although here we hnd appreciable cooling only close to the Eddington 
rate, we expect that for near-Keplerian accretion hows Compton scattering will become dy¬ 
namically important at lower accretion rates, as individual huid elements will have more 
time to cool (their radial velocity is lower) before accreting onto the black hole. 


5.2. Axisymmetric Radiating Kerr Black Hole Accretion 

As another preliminary application of bhlight we consider the ehect of radiation on 
an intermediate accretion rate black hole accretion how. Recall that for systems with 
L ~ 10“®LEdd (like Sgr A*) radiation will have little ehect since the cooling timescale is 
long compared to the accretion timescale. This is the classical RIAF regime. For systems 
with L ~ 10“®LEdd (like M87) the ehect of radiation depends on temperature and the dis¬ 
tribution function of the electrons, but in certain circumstances radiation interactions - 
especially Compton cooling - can cool the how on timescales comparable to or shorter than 
the accretion timescale. 

Table 4. Parameters for radiating Bondi accretion 


Case 

n 

cm“^ 

m 

L 

L'Edd 

1 

3.0 X 10® 

4.01 X 10-^ 

2.03 X 10-1^ 

2 

3.0 X 10^ 

4.01 X 10“® 

1.46 X 10-12 

3 

3.0 X 10® 

4.01 X 10“® 

2.08 X 10-1® 

4 

3.0 X 10*^ 

4.01 X 10“^ 

5.81 X 10“^ 

5 

3.0 X 10^° 

4.01 X 10“® 

3.05 X 10-^ 
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Fig. 16.— Fluid and radiation profiles for the radiating Bondi problem. Case 1 is shown in 
purple, Case 2 in teal, Case 3 in red. Case 4 in green, and Case 5 in blue. The dashed line 
shows Case 5 without Compton scattering. 
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Fig. 18.— Self-convergence for the radiating Bondi problem near the Eddington limit. 
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Here we consider an axisymmetric accretion flow with parameters inspired by M87. We 
set black hole mass M = 6.6 x IO^Mq following Gebhardt et al. (2011) and dimensionless 
spin a* = 0.9375, and we adjust Ai so that m ~ 6.3 x 10“® (Kuo et al. 2014). Although the 
distribution functions for ions and electrons in M87-like systems probably contains multiple 
components, we adopt thermal distributions for both and set the ion and electron tempera¬ 
tures Ti = STg. Ti/Tf. may strongly affect the dynamics of this system; essentially, it controls 
the cooling rate. We set the adiabatic index 7 = 13/9, appropriate for ionized hydrogen 
when kTi rripC^ and kTe 3> rrieC^- We include synchrotron emission for a relativistic, 
thermal distribution of electrons (emissivity given by Equation 72 of Leung et al. 2011), 
thermal absorption, and Compton scattering. Bremsstrahlung is neglected, as it is a small 
correction to the emissivity within ~ lO^M of the black hole for such intermediate accretion 
rate systems (e.g. Narayan & Yi 1995). 

The initial conditions are a Fishbone-Moncrief torus (Fishbone & Moncrief 1976) with an 
inner radius at QGM/c^ and pressure maximum at 12GM/c^. We extend this configuration 
by adding a weak poloidal magnetic held that follows isodensity contours, using the same 
procedure as in McKinney & Gammie (2004), but with the vector potential modified by a 
cos 6 factor to produce a two-loop configuration. The held strength is normalized so that (3 
has a global minimum value of 100. Small perturbations are applied to the internal energy 
to efficiently initiate the magnetorotational instability. No radiation is present in the initial 
conditions. 

We use the modified Kerr-Schild coordinates (see §3.5), with ro = 0 and hg = 0.3. The 
grid runs from r = 0.98(1 -|- a/1 — al)GM/c^ to r = 40GM/c^, and from 6 ^ = 0 to 6 * = tt. We 
impose outflow boundary conditions on both the fluid and radiation at the inner and outer 
radial boundaries, and reflecting polar boundary conditions. In this instance we adopt the 
piecewise parabolic method for reconstruction. We set the GFL number to 0.7, and evolve 
the system until t = 2000 M. 

To accurately sample cooling due to scattering events, we bias the scattering via the 
method given in Section 3.2.7 such that hg oc 0g. To determine the absolute number of 
superphotons in steady state, we require only that the bolometric light curve be satisfactorily 
resolved (here, ~ 1.1 x 10® superphotons at any one time). In the future we will more carefully 
consider the resolution requirement for these models. 

We find qualitative differences between the radiative torus simulated in bhlight and 
the same model evolved with ideal GRMHD due mainly to Gompton cooling in the hot, 
dense regions of the flow. Figure 19 shows fluid temperature and comoving radiation energy 
density at t = 2000M. Figure 20 shows density contours at t = 2000M for the two models. 
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Figure 21 shows shell-averaged density-weighted temperature, 

/o \ - / 

' “ I V^dxWp ’ 


(61) 


as a function of radius for both bhlight and ideal GRMHD calculations; evidently the 
plasma is signihcantly cooler in the model with radiation interactions. Notice, however, that 
onr single temperatnre model may be having an nnrealistically strong dynamical effect; at 
these accretion rates the ions are likely imperfectly conpled to the electrons. 


We evalnate the Inminosity L, Eqnation 59, at large radius by integrating over all zones 
in a spherical shell. Figure 22 shows L and the radiative efficiency r/ = L/Mc^. The mean 
T] we hnd between t = 1750M and t = 2000M, ( 77 ) = 0.51, where 


/ \ 

{v) = 


f dt Mc^' 


(62) 


is high compared to the thin disk value, ~ 0.18 (for this a*). Most of this energy is extracted 
by Compton scattering at r ~ 10 — 15GM/c^. This high efficiency is a consequence of the 
flow having not yet reached steady state. 

We plan to explore radiative flows at intermediate accretion rates with more sophisti¬ 
cated treatments of the electron physics in fntnre work. 
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Fig. 19.— Torus temperature ©e and comoving radiation energy density at t = 

1500M. 
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Fig. 20.— Comparison of gas density p between GRMHD and bhlight torus calculations. 
Note especially that for the bhlight result, the disk is relatively thin and the funnel is poorly 
developed. 
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Fig. 22.— Torus luminosity, instantaneous efficiency 77 , and mass accretion rate as a function 
of time. The dashed line denotes the thin disk efficiency at this spin. The grey region 
indicates the portion of rj prior to the onset of accretion across the horizon which we omit. 
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6. Conclusion 

We have introduced bhlight which, in coupling a second order Godunov scheme to 
a frequency-dependent Monte Carlo radiative transfer scheme, provides a full solution to 
the equations of GRRMHD. bhlight displays convergence on a number of test problems, 
and we have demonstrated evolution of our target application: relativistic accretion flows at 
moderate accretion rate, bhlight enables more nearly ab initio study of flows in this regime, 
which have historically proven resistant to other methods of inquiry. 

Numerical schemes have limitations. Apart from failure modes related to large optical 
depths and large radiation pressures mentioned previously (which prohibit near-Eddington 
studies as of this writing), Monte Carlo techniques are simply expensive, particularly when 
geodesics are nontrivial. The axisymmetric torus problem we reported was performed on 
one 8-core node (using mpirun to alternate between fluid and radiation MPI sectors) for 146 
hours, achieving approximately 9.5 x 10"^ superphoton updates (interactions and geodesic 
steps) per core-second. In comparison, the pure GRMHD torus run required only 69 core¬ 
hours, about 17 times less expensive even at this low superphoton resolution. The true 
minimum relative cost of bhlight over harm, however, will depend on the required resolution 
in the specific intensity for the particular problem at hand. 
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A. Radiation Boundary Conditions 

In the general coordinate frame, certain boundary conditions on the superphotons re¬ 
quired for test problems in bhlight are more complex than in the nonrelativistic case. Here 
we review such boundaries as used. 


A.l. Reflecting Boundary Conditions 

The static atmosphere test in Section 4.7 uses reflecting boundaries that are aligned with 
the coordinates. How should one apply the reflecting boundary conditions to the radiation 
held? This is not trivial in Kerr-Schild coordinates where the naive approach of simply 
changing the sign of the radial component of the wavevector is wrong, because the radial 
component of the shift does not vanish. 

When a photon crosses the boundary, we build an orthonormal tetrad with time com¬ 
ponent = (m*, 0, 0 , 0) (i.e. the tetrad is stationary in the coordinate, and hence the 
boundary, frame) and one spatial component that is normal to the boundary. We transform 
the wavevector to the tetrad frame, reverse the normal component of the wavevector, and 
transform back to the coordinate frame. 


A.2. Equilibrium Boundary Conditions 

For problems with fluid inflow across the boundary (e.g. the relativistic shocks in Section 
4.6) the fluid advects a thermal radiation held with it across the boundary. How should one 
sample the incoming photons on the boundary? The problem is that one is sampling a flux 
rather than the distribution function itself. 

We have found that the simplest procedure is to sample the distribution function and 
multiply the weights in the sample by cos 6*, where 6 is the angle between the wavevector 
and the normal to the boundary. 


B. The Radiating Atmosphere Under Full Transfer 

We revisit the analytic solution to the radiating atmosphere problem described in Sec¬ 
tion 4.7, in which fluid and radiation conhned in the Schwarzschild spacetime between two 
reflective spherical shells maintain a static atmosphere, without resorting to any closure for 
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describing the radiation. 

Consider the relativistic transfer equation in invariant form (ignoring scattering), 


D 




(Bl) 


which may be rewritten in terms of the invariant differential optical depth dr = ua^^dX (and 
noting that the Planck function Bi, = ipjay, i.e. Kirchoff’s law applies) to give 


A (hd 

dr V z/3 


5 , 






(B2) 


where all quantities in brackets are again invariant. 

We adopt the ansatz for the fluid temperature distribution y/—gooT{r) = Too from 
Section 4.7. Consider also the frequency u = —k^u^ of a superphoton at radius r in the 
Schwarzschild spacetime. The four-velocity of a local frame not moving with respect to the 
coordinate system is = (1/^^— 5 ^ 00 , 0, 0, 0): then, u = —/co/V—fi'oo- The invariant Planck 
function /c^) f {In'/ksT). However, hv/ksT = —/i/co/Stt/cbToo is constant along 

geodesics; Byj is therefore also constant along every ray. 

For reflecting boundary conditions, rays of the intensity 1^ do not terminate; they 
instead repeatedly reflect off the boundaries all the way back to dr — )■ 00. Because By/v^ is 
constant along every ray, for the stationary system we have ly/u^ = By/v^ everywhere. Our 
assumed temperature distribution is therefore consistent with the full transfer equation, and 
the solution presented in Section 4.7 is exact for all optical depths. 


C. Linear Modes in Relativistic Radiation Magnetohydrodynamics 

We present the linearized equations of radiation magnetohydrodynamics in flat space- 
time, assuming the Eddington closure (Farris et ah 2008) with a gray absorption opacity k 
and setting ks = c = 1. The governing equations are then given in terms of divergences of 
the matter four-current and the MHD and radiation stress-energy tensors, along with the 
magnetic induction equation. In plane-parallel symmetry, we search for wave solutions of 
the form P = (p, m, v}, v?, B^, E, T^) = (po -|- 6p, uq -f 6u, 6u^, Bq, Bq + 6B‘^, Eq + 
SE, SF^, 5F‘^) (i.e. we conhne variation to a plane), where 5 oc exp(ci;t — ikx) is a small per¬ 
turbation, and Eq = a/j ((7 — 1)mo/po)"^ to enforce radiative equilibrium of the background 
state. We write the linearized systems in the form cahP = A(5P; the dispersion relation is 
then det(A — lea) = 0. We hnd that the matrix A is 
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AikBlEQ{'y-l) 

0 

0 

—AikBQ Eq 
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-iKpoB^Eo 

—C 3 KP 0 

3Ci 

3(C2+B2) 

3Ci 
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where 

C\ = 7^Mo + 2 ( 7 % + Po)Bq + 2'ypQUo + pg, 

C 2 = Bq + 7M0 + po, 

(73 = aCi + 4E0C2. 
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